#include "Geometry.h"
using namespace blitz;

// Open Geometry File
void TriangleMesh::ReadGeometry(char* GeoName,emx::Real3D Shift)
{
	std::cout<<"Opening Geometry \""<<GeoName<<"\""<<endl;
	fin.open(GeoName);									
	if (!fin) std::cerr<<"File not exist!"<<std::endl;
	ReadNodes(Shift);
	ReadPatchs();
	ReadEdges();
	fin.close();

	//Saving Edge Length
	for (int i=0; i<Edge_Num; i++)
	{
		emx::Real3D vec = Edge(i)->EndVertex->Position-Edge(i)->BeginVertex->Position;
		Edge(i)->EdgeLength = sqrt(dot(vec,vec));
	}

	//Saving Patch Center, Area and Normal Vector
	for (int i=0; i<Patch_Num; i++)
	{
		int node1 = Patch(i)->Vertex(0)->Idx;
		int node2 = Patch(i)->Vertex(1)->Idx;
		int node3 = Patch(i)->Vertex(2)->Idx;
		emx::Real3D R1,R2,R3;
		R1 = Node(node1)->Position;
		R2 = Node(node2)->Position;
		R3 = Node(node3)->Position;
		Patch(i)->Center = (R1+R2+R3)/3;
		Patch(i)->Point(0) = Patch(i)->Center;
		Patch(i)->Point(1) = R1/6.0+R2/6.0+2.0*R3/3.0;
		Patch(i)->Point(2) = 2.0*R1/3.0+R2/6.0+R3/6.0;
		Patch(i)->Point(3) = R1/6.0+2.0*R2/3.0+R3/6.0;
		emx::Real3D R12,R23,R31;
		R12 = R2-R1;
		R23 = R3-R2;
		R31 = R1-R3;
		emx::Real3D tmp = cross(R12,R23);
		Patch(i)->Normal = tmp/sqrt(dot(tmp,tmp));
		emx::real_data e1 = sqrt(dot(R12,R12));
		emx::real_data e2 = sqrt(dot(R23,R23));
		emx::real_data e3 = sqrt(dot(R31,R31));
		emx::real_data s = 0.5*(e1+e2+e3);
		Patch(i)->Area = sqrt(abs(s*(s-e1)*(s-e2)*(s-e3)));
	}

	//Saving Rho^plus and Rho^minus
	for (int i=0; i<Edge_Num; i++)
	{
		//Saving vector \rho^plus
		int node3 = Edge(i)->LeftFacet->Vertex(0)->Idx+Edge(i)->LeftFacet->Vertex(1)->Idx+Edge(i)->LeftFacet->Vertex(2)->Idx
			-Edge(i)->BeginVertex->Idx-Edge(i)->EndVertex->Idx;
		Edge(i)->LeftPoint = Node(node3)->Position;
		Edge(i)->Rho_plus = Edge(i)->LeftFacet->Center - Node(node3)->Position;
		//Saving vector \rho^minus
		int node4 = Edge(i)->RightFacet->Vertex(0)->Idx+Edge(i)->RightFacet->Vertex(1)->Idx+Edge(i)->RightFacet->Vertex(2)->Idx
			-Edge(i)->BeginVertex->Idx-Edge(i)->EndVertex->Idx;
		Edge(i)->RightPoint = Node(node4)->Position;
		Edge(i)->Rho_minus = Node(node4)->Position - Edge(i)->RightFacet->Center;
	}
}

// Import Node Information
bool TriangleMesh::ReadNodes(emx::Real3D Shift)										
{
	fin>>Node_Num;															
	Node.resize(Node_Num);
	int idx;
	for (int i=0; i<Node_Num; i++)
	{
		Node(i) = new node;
		Node(i)->Idx = i;
		fin>>idx>>Node(i)->Position[0]>>Node(i)->Position[1]>>Node(i)->Position[2];
		Node(i)->Position += Shift;
	}
	
	return true;
}

// Import Patch Information
bool TriangleMesh::ReadPatchs()									
{
	fin>>Patch_Num;														
	Patch.resize(Patch_Num);
	int idx, node1, node2, node3;
	for (int i=0; i<Patch_Num; i++)
	{
		Patch(i) = new facet;
		Patch(i)->Idx = i;
		fin>>idx>>node1>>node2>>node3;
		Patch(i)->Vertex(0) = Node(node1-1);
		Patch(i)->Vertex(1) = Node(node2-1);
		Patch(i)->Vertex(2) = Node(node3-1);
	}

	return true;
}

// Import Edge Information
bool TriangleMesh::ReadEdges()										
{
	fin>>Edge_Num;
	Edge.resize(Edge_Num);
	int idx, node1, node2, leftpat, rightpat;
	for (int i=0; i<Edge_Num; i++)
	{
		Edge(i) = new edge;
		Edge(i)->Idx = i;
		fin>>idx>>node1>>node2>>leftpat>>rightpat;
		Edge(i)->BeginVertex = Node(node1-1);
		Edge(i)->EndVertex = Node(node2-1);
		Edge(i)->MidPoint = (Edge(i)->BeginVertex->Position+Edge(i)->EndVertex->Position)/2;
		Edge(i)->LeftFacet = Patch(leftpat-1);
		Edge(i)->RightFacet = Patch(rightpat-1);
	}

	return true;
}
